home *** CD-ROM | disk | FTP | other *** search
/ Best of www.BestZips.com (Collector's Edition) / Best of WWW.BESTZIPS.COM Collector's Edition (JCSM Shareware) (JCS Marketing).ISO / prgtools / tn2.zip / LEAST_SQ.T < prev    next >
Text File  |  1996-11-15  |  1KB  |  77 lines
  1. %
  2. % "least_sq.t" solves for least square error
  3. % coefficients in:
  4. %
  5. %       y = a + b x + c x^2
  6. %
  7. %   Sample program for the T Interpreter by:
  8. %
  9. %   Stephen R. Schmitt
  10. %   962 Depot Road
  11. %   Boxborough, MA 01719
  12. %
  13.  
  14. const DIM : int := 3
  15.  
  16. program
  17.  
  18.     var n : int
  19.     var x, y , det: real
  20.     var v, w : rmatrix
  21.     var a, b : rvector
  22.     label program_exit :
  23.  
  24.     prompt "how many data points?"
  25.     get n
  26.  
  27.     put "           x           y"
  28.     loop
  29.  
  30.         exit when n <= 0
  31.  
  32.         prompt "x?"
  33.         get x
  34.         prompt "y?"
  35.         get y
  36.  
  37.         put x:12:6, y:12:6
  38.  
  39.         w[0,0] := w[0,0] + 1.0
  40.         w[0,1] := w[0,1] + x
  41.         w[0,2] := w[0,2] + x * x
  42.  
  43.         w[1,0] := w[1,0] + x
  44.         w[1,1] := w[1,1] + x * x
  45.         w[1,2] := w[1,2] + x * x * x
  46.  
  47.         w[2,0] := w[2,0] + x * x
  48.         w[2,1] := w[2,1] + x * x * x
  49.         w[2,2] := w[2,2] + x * x * x * x
  50.  
  51.         b[0] := b[0] + y
  52.         b[1] := b[1] + x * y
  53.         b[2] := b[2] + x * x * y
  54.     
  55.         decr n
  56.  
  57.     end loop
  58.  
  59.     det := invert( w, v, true )
  60.  
  61.     if det = 0.0 then
  62.  
  63.         put "singular!"
  64.         goto program_exit
  65.  
  66.     end if
  67.  
  68.     mul_mat_vec( v, b, a )
  69.  
  70.     put "\ncoefficients:"
  71.     put "a = ", a[0]:12:8
  72.     put "b = ", a[1]:12:8
  73.     put "c = ", a[2]:12:8
  74.  
  75.     program_exit:
  76.  
  77. end program